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Online  vehicle  trajectory  reshaping  is  desired  for  a  class  of  autonomous  air 
vehicles  such  as  RLVs  in  order  to  avoid  catastrophic  failure  when  subjected  to 
performance  restricting  damages  and  failures.  An  Adaptive  Trajectory  Reshaping 
and  Control 1  (ATRC)  system  is  envisioned  that  responds  to  altered  vehicle  conditions 
by  continuously  retargeting  and  reshaping  the  reference  RLV  trajectory  satisfying 
the  feasibility  constraints.  On-line  trajectory  reshaping  to  determine  a  feasible 
reference  trajectory  is  computationally  a  difficult  problem  for  real  time 
applications.  ATRC  is  exploring  the  principles  of  vehicle  dynamics  inversion  for  on¬ 
line  generation  of  feasible  reference  trajectory.  Two  essential  components  for 
generating  reference  trajectory  for  air-vehicles  using  “inverse  dynamics” 
methodology  are  aerodynamic  model  of  the  vehicle  that  is  representative  of  the 
current  state  of  the  vehicle,  and  a  framework  for  modeling  the  vehicle  trajectory. 
Physics  based  modeling  such  as  DATCOM  allows  fast  computation  of  aerodynamic 
coefficients  for  given  flight  points  and  the  results  can  be  stored  in  a  tabular  form. 
However,  for  efficient  real-time  trajectory  reshaping  application,  it  is  desired  to 
represent  aerodynamic  coefficients  in  smooth  functional  form  that  are  governed  by 
few  parameters.  Similarly,  trajectories  must  also  be  represented  by  smooth 
functions.  In  this  paper  we  present  modeling  of  smooth  functions  using  a  set  of 
basis  functions  that  are  suitable  for  trajectory  reshaping  of  the  air  vehicles.  A 
desirable  feature  for  function  modeling  is  the  easy  imposition  of  boundary  as  well  as 
mid  point  constraints  in  the  function  using  small  number  of  parameters  without 
limiting  the  scope  of  the  function.  In  this  paper  we  present  a  design  of  polynomial 
based  set  of  constrained  orthonormal  basis  functions  in  one  and  two  dimension. 


I.  Introduction 

The  large  potential  for  space  utilization  is  not  being  exploited  as  it  is  currently  inhibited  by  the  huge 
cost  of  launching  operations.  The  benefit  of  advance  space  utilization  can  be  greatly  increased  by  making 
space  utilization  more  affordable  The  Reusable  Launch  Vehicle  (RLV)  programs  are  targeted  towards 
affordable  space  utilization.  However,  to  maintain  the  economical  viability  of  RLVs,  it  is  important  to 
enhance  operations  safety  and  reliability  by  providing  the  RLV  the  capability  to  respond  to  various 
uncertainties  and  emerging  emergency  situations.  Responding  to  an  uncertain  environment  after  a 
damage/failure  presents  many  tough  technical  problems  for  this  class  of  vehicle.  These  problem  manifests 
in  the  following  challenges  that  must  be  overcome:  First,  to  adequately  determine  and  model  the  dynamic 
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characteristics  of  the  vehicle  in  the  altered  state  after  a  damage/failure;  second,  to  estimate  the  new 
constraints  and  limitation(s)  of  the  vehicle;  third,  to  adopt  and  reconfigure  the  command,  control  and 
guidance  of  the  vehicle  to  the  modified  system  dynamics;  and  fourth,  to  design  and  plan  a  new  feasible  path 
with  respect  to  the  end  goal  maximization.  A  high  percentage  of  such  damage/failure  cases  leave  the 
vehicle  in  an  uncontrolled  and  uncertain  environment  with  a  highly  likely  probability  of  ultimately  entering 
a  state  for  catastrophic  failure.  The  high  cost  of  loss  resulting  from  catastrophic  failures,  has  prompted 
researchers  in  the  direction  of  developing  technologies  to  assist  in  minimizing  such  failures. 

Damage  to  a  vehicle  or  a  sub-system  failure  may  result  in  modification  of  the  applicable  trajectory 
constraints  and/or  the  dynamical  behavior  of  the  system,  consequently  making  the  previously  designed 
reference  trajectory  infeasible.  An  acceptable  trajectory  for  a  dynamical  system  is  a  solution  of  a  two-point 
boundary  value  problem  for  a  set  of  governing  differential  equation  of  motion.  The  real  world  systems  such 
as  air  vehicles  are  highly  non-linear  systems  and  impose  a  set  of  constraints  on  the  trajectory  variables  as 
well  as  control  variables.  In  inverse  dynamics  approach,  a  trajectory  is  specified  first,  which  results  in 
solving  a  set  of  algebraic  equations,  yet  strictly  satisfying  the  non-linear  differential  equations  of  a  non-flat 
system. 

In  this  paper  first  we  introduce  the  architecture  of  an  Adaptive  Trajectory  Reshaping  and  Control 
(ATRC)  system  for  general  class  of  RLV  system,  which  is  based  on  the  principles  as  described  in  Ref.  1 . 
Next  we  discuss  the  general  inverse  dynamics  approach  for  trajectory  reshaping  of  RLVs  and  the 
motivation  for  functional  modeling.  For  real  time  trajectory  determination,  there  are  two  types  of  functions 
that  must  be  modeled.  One  set  of  function  is  needed  to  describe  the  vehicle  trajectory.  The  functions  for  the 
spatial  coordinates  of  the  vehicle  constitute  a  vehicle  trajectory.  These  trajectory  functions  are  normally 
independent  in  one  variable.  The  second  set  of  functions  that  must  be  modeled  in  real  time  is  to 
approximate  the  modified  aerodynamic  constraints.  The  aerodynamic  constraints  are  modeled  by  defining 
aerodynamic  coefficients  for  the  flight  envelope  of  the  vehicle.  Mostly,  these  functions  are  two  dimensional 
with  Mach  and  angle  of  attack  being  the  independent  parameters.  A  general  approach  to  model  a  function 
is  to  parameterize  the  function  and  then  determine  the  parameters  that  satisfy  any  governing  constraints  to  a 
satisfactory  level.  In  section  IV,  we  describe  a  novel  approach  to  design  a  set  of  basis  functions  with  a  class 
of  constraints  built  in,  that  reduces  the  complexity  of  the  solution  and  results  in  better  approximations. 


II.  Adaptive  Trajectory  Reshaping  and  Control  (ATRC)  System 

The  ATRC  system  enhances  RLV  capability  to  avoid  catastrophic  failure  when  subjected  to 
performance  restricting  damages  and  failures.  The  overall  goal  of  ATRC  translates  into  specific 
requirements  for  design  and  development  of  functionalities  related  to  adaptable  and  reconfigurable 
command,  control,  and  guidance  system  for  the  RLVs.  and  real  time  solution  techniques  for  an. 

Figure  1  shows  the  general  architecture  of  the  envisioned  ATRC  system  for  RLVs.  Note  that  the  above 
structure  is  specific  to  longitudinal  motion  of  the  vehicle,  however,  it  can  be  easily  extended  to  six  degree 
of  freedom.  The  main  components  of  the  envisioned  ATRC  system  requires: 

1 .  On-line  system  identification  that  includes  parameter  estimation  and  parameter  projection  for 
constraint  boundary  determination.  The  constraint  boundaries  influence  the  trajectory  reshaping  of  the 
vehicle. 

2.  Real  time  trajectory  determination  for  reshaping  reference  trajectory  under  feasibility  constraints. 

3.  Adaptive,  closed  loop  control  and  guidance  system  for  reference  trajectory  tracking. 
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III.  Inverse  Dynamics  Approach  for  Trajectory  Generation 

The  advantage  of  inverse  dynamics  approach  is  that  one  can  avoid  integration  of  differential  equations 
altogether.  The  main  difference  is  that  we  parameterize  the  trajectory  itself  and  then  uses  numerical 
techniques  for  solving  these  parameters  that  minimizes  the  objective  function  and  other  inequality 
constraints.  Once  a  smooth  trajectory  is  specified,  the  time  derivatives  are  also  fixed  and  hence  the 
differential  equations  become  simple  algebraic  equations. 

A.  Reference  Trajectory  Design 

To  determine  a  trajectory  for  a  non-linear  dynamic  system,  a  solution  must  be  found  that  satisfies  the 
set  of  differential  equations  governing  the  dynamics  of  the  system.  Further,  the  trajectory  solution  should 
not  violate  some  non-linear  constraints,  which  limits  the  operation  capability  of  the  system.  For  an  aircraft, 
the  constraints  arise  due  to  several  factors  such  as  the  limitations  on  angle  of  attack ,  load  factor , 
aerodynamic  heating ,  and  actuator  saturation. 

There  are  two  primary  approaches  for  trajectory  generation  and  these  have  been  classified  in  the 
literature  [8]  as  the  “integral  approach”  and  the  “differential  approach.”  In  any  approach,  where  generation 
of  a  trajectory  involves  the  integration  of  the  equations  of  motion,  this  approach  is  classified  as  the 
“integral  approach.”  In  a  differential  approach,  an  assumed  functional  form  for  trajectory  is  differentiated 
to  obtain  algebraic  functions  for  the  higher  derivatives,  which  are  required  to  impose  constraints  on  the 
control  inputs  for  the  “inverse  dynamics”  solution.  There  are  various  applications  where  inverse  dynamics 
have  been  used,  such  as  spacecraft  trajectories  and  path  planning  in  robotics  [9]  and  overhead  cranes  [10]. 
Historically,  the  inverse  dynamics  approach  has  been  used  for  “differentially  flat”  systems.  A  system  is 
“differentially  flat”  [see  11,  12]  if  there  exists  a  set  of  outputs,  known  as  “flat  outputs,”  such  that  there  is  a 
one-to-one  correspondence  between  the  trajectories  of  flat  outputs  and  the  full  state  and  control  inputs  of 
the  system.  In  our  approach,  we  use  an  inverse  dynamics  approach  for  determining  trajectory  as  this  allows 
us  to  solve  algebraic  equations  instead  of  integrating  ODEs.  With  the  inverse  dynamics  approach  for 
aircraft  trajectories,  a  problem  arises  due  to  inherent  under- actuation  in  most  of  the  aircrafts.  For  a  six- 
degree  of  freedom  aircraft,  there  are  normally  four  controls:  thrust,  elevator,  aileron,  and  rudder.  Ref.  [7] 
defines  a  novel  trajectory  generation  scheme,  which  uses  pseudo  forces  for  inverse  dynamic  computation. 
In  this  paper  we  concentrate  on  function  modeling  approach  that  to  make  the  inverse  dynamic  solution 
approach  more  efficient. 
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B.  Inverse  Dynamics  Approach 

Assume  X  be  the  state  vector  for  a  vehicle  and  X  =  x(Y)  represent  the  trajectory  of  the  vehicle.  The 
governing  equations  of  motion  for  an  air  vehicle  can  be  represented  in  the  following  form 

Ml 

x  =  \C (v  (x))]<  x  >  =  Cj  (v(x))  +  C2  (v(x))x  +  C3  (v  (x))i/ ,  ( 1 ) 

u 

where  v(x)  are  the  independent  parameters  defined  in  terms  of  the  subset  of  vehicle  states,  and 
Cf(v(x))  are  the  non-linear  coefficient  functions  in  independent  parameters. 

For  example,  the  pitch  dynamics  of  a  vehicle  for  longitudinal  dynamics  can  be  written  as 

lq  =  CnjM,a)+Cnlj  (M ,a)q  +  Cm&  (M ,a)Se  .  (2) 

The  inverse  dynamic  solution  for  Eq.  (1)  can  be  written  as 

U  =  C3_1(v(x)Xx-C1(v(x))+C2(v(x))x)  (3) 

For  a  real  time  solution,  in  presence  of  uncertainty  due  to  damage  to  the  vehicle  requires  modeling  of 
the  trajectory  x  =  x(r)  and  estimated  models  of  the  coefficient  functions  Cz  (v(x)) .  In  the  following  section 
we  will  address  one  dimensional  trajectory  modeling  and  two-dimensional  coefficient  functional 
representation  in  parametric  form. 

IV.  One-Dimensional  Function  Modeling  with  Constraints  using  Basis  Functions 

Often  we  need  to  estimate  or  model  non-linear  functions  that  must  satisfy  some  constraints.  For 
example,  the  inverse  dynamic  solution  approach  requires  parameterization  of  multiple  system  outputs 
variables  or  the  trajectory  variables  with  some  constraints  imposed  on  it.  Note  that  the  feasible  trajectory 
must  also  satisfy  the  governing  differential  equations.  The  trajectory  parameters  are  determined  using  an 
optimization  process  that  minimizes  the  error  between  trajectory  variables  and  the  solution  of  the  governing 
equations.  Note  that  if  some  of  the  trajectory  constraints,  such  as  the  boundary  constraints  are  forced  in  the 
trajectory  parameterization  itself,  it  helps  in  reducing  the  complexity  of  the  optimization  problem.  Verma  et 
al  [4,7]  presented  the  following  approach  to  represent  trajectory  that  is  useful  to  impose  boundary  and 
inpoint  constraints.  To  prescribe  a  smooth  trajectory  functional  for  a  specific  trajectory  coordinate,  the 
position  coordinate  is  represented  by  a  twice  differentiable,  smooth  function  so  that  velocity  and 
acceleration  can  be  uniquely  defined.  The  path  for  an  individual  position  coordinate  is  defined  as  a 
function  of  normalized  time  r-t/Tf ,  where  Tf  is  the  total  time  for  the  maneuver.  An  individual 
coordinate  trajectory  is  structured  to  have  two  parts,  a  base  trajectory  and  a  perturbation  of  the  base 
trajectory.  If  P(r)  is  the  ith  coordinate,  it  is  chosen  to  be  of  the  form 

P{t)  =  P0  (r)  +  P{  {tTLUjVj  (t)  (4) 

7=1 

Here  P0(r)  is  the  base  trajectory  function,  which  is  chosen  such  that  it  satisfies  any  boundary  or  mid 
point  constraints  for  the  the  trajectory.  One  example  approach  P0(r)  could  be  chosen  as  a  minimum  degree 
polynomial  spline  that  satisfies  the  desired  constraint  conditions.  The  base  trajectory  by  itself  may  not  be  a 
true  solution  of  the  governing  equations  and  hence  an  infeasible  trajectory.  To  make  the  trajectory  feasible, 
we  add  a  second  term  on  the  right  hand  side  of  the  trajectory  equation.  The  perturbation  term  must  be 
designed  to  ensure  that  the  overall  trajectory  function  P{r)  is  the  solution  of  the  governing  equations. 
Perturbation  term  is  defined  using  a  set  of  ortho-normal  basis  functions  <Pj(t)  that  are  modified  by  the 
term  Px{r) .  The  term  Px{r)  is  a  weight  polynomial  constraining  the  perturbation  term  to  contribute  zero 
value  for  the  already  satisfied  desired  constraints  by  base  trajectory  P${r).  For  example  P\{r)  can  be 
chosen  as 

Pi(r)  =  (r-a)/,(T-b)g ,  (5) 

where  p  and  q  are  integers,  usually  less  or  equal  to  3.  The  integers  p  and  q  depend  upon  the  highest 
order  condition  match  required  at  times  r  =  a  and  r  =  b  respectively.  For  example,  when  boundary 
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conditions  at  r  =  0  and  r  =  1  up  to  acceleration  level  must  be  imposed  on  the  trajectory  function,  Px{r)  can 
be  chosen  as  r3(r-l)3 .  The  drawback  of  this  approach  is  that  the  after  basis  functions  (pj  are  modified  by 

Px  (r) ,  it  looses  the  orthoganality,  which  makes  it  computationally  harder  to  determine  the  coefficients  for 
the  trajectory  solution.  In  the  next  section  we  present  an  approach  to  design  orthogonal  polynomial  basis 
functions  with  inherent  desired  constraints. 

V.  Constraint  Orthogonal  Polynomial  Basis  Functions 

The  choice  of  basis  functions  ultimately  influences  the  complexity  of  the  modeling  of  a  desired  function 
or  behavior.  It  is  well  known  that  a  sub  set  of  a  complete  set  of  basis  functions  can  approximate  any  given 
function  to  a  desired  accuracy  by  choosing  a  sufficient  number  of  basis  functions  elements  in  the  sub  set. 
However,  if  a  function  to  be  approximated  must  meet  certain  constraints,  it  presents  various  problems. 
First,  a  finite  set  of  basis  function  may  not  ensure  the  constraints  satisfaction  on  the  function  to  be 
approximated.  Second,  a  large  number  of  basis  functions  may  be  required  for  the  satisfactory  constraint 
approximation.  Now,  if  we  can  create  a  set  of  basis  functions,  where  all  the  basis  functions  satisfy  the  given 
constraints,  any  linear  combination  of  those  functions  will  also  satisfy  the  given  constraints.  With  this 
motivation  we  developed  an  approach  to  design  a  constraint  orthogonal  polynomial  basis  functions.  To 
explain  our  approach,  we  first  present  a  way  for  generating  one-dimensional  unconstraint  orthogonal 
polynomial  based  basis  functions. 

A.  One  Dimensional  Basis  Functions 

Let  ®,(x)  represent  the  1-D  zth  function  of  the  orthonormal  polynomial  basis  functions  given  as: 

<$>i(x)=Yjakxk,  i  =  OX-  (6) 

k=0,i 

Notice  that  the  zth  basis  function  has  /+ 1  polynomial  coefficients  to  be  determined.  We  impose  the 
constraints  on  the  functions  to  determine  the  coefficients  to  ensure  orthogonality  and  normalization. 
Defining  the  inner  product  of  two  functional  as 

l 

(<!>,■  (x),<D*(x))  =  J  ®,(x)c Pk(x)dx .  (7) 

0 

The  norm  and  the  orthogonality  conditions  can  be  obtained  as: 

l 

J  (O;  (x))2dx  =  1 ,  Normalizing  condition.  (8) 

o 
i 

Joz(x)O^(x)ix  =  0,  k  =  Orthogonality  condition.  (9) 

o 

The  orthogonality  condition  gives  i  constraints  while  an  additional  constraint  from  normalization  ensures 
that  all  the  coefficients  can  be  obtained. 

For  automatically  generating  the  basis  set  we  used  the  following  recursive  approach.  Assuming  that  we 
already  know  0>k,k  =  0,/ -1 .  Define  an  arbitrary  polynomial  /J(x)  of  degree  i.  For  k  =  0  to  i  -  1  we 

modify  fi(x)  as 

/*+1(x)  =  /*(x)-(  ®*(x)/*(x))®t(x)  (10) 

The  ith  basis  function  ^>t(x)  is  now  given  as 

*/(*)  =  //(*)  (11) 

For  constrained  orthogonal  polynomial  function,  Eq.  (  6)  is  modified  as  follows. 

(12) 

k=0,i 
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where,  y/{x)  imposes  a  desired  constraint  on  the  basis  function.  Figure  2  and  Figure  3  show  some 
examples  of  a  class  of  orthogonal  basis  functions  that  incorporate  the  given  constraints.  The  example 
presents  normalized  orthogonal  basis  functions  in  the  range  0  to  1 .  In 


x  - ►  -  x  - ► 


(a)  (b) 

Figure  2.  Constraint  Polynomial  Basis  Functions,  (a)  Basis  Function  are  Constraint  to  be  Zero  atx  = 
0.5.  (b)  Both,  Function  and  its  Derivative  are  Constrained  to  be  Zeros  atx  =  0.5 

Figure  2(a),  the  value  of  the  function  is  constraint  to  be  zero  mid  way,  i.e.  /(x)|x=05  =  0  Figure  2(b)  adds 
an  additional  constraint  on  the  slope  of  the  curve  or  the  first  derivative  of  the  function  given  as 
df(x)/dx\x=05  =  0  .  Note  that  the  first  set  of  constraint  results  in  anti-symmetric  orthogonal  basis  functions 

as  can  be  seen  in  Figure  2(a),  while  the  second  set  of  constraint  results  in  symmetric  orthogonal  basis 
functions  as  shown  in  Figure  2(b).  In  Figure  3,  the  application  of  the  constraint  has  been  shifted  to  the 
quarter  point,  i.e.  x  =  0.25  .  In  this  case,  the  basic  nature  of  orthogonal  basis  functions  remains  the  same, 
albeit  with  a  loss  of  symmetry.  Figure  2  and  Figure  3  demonstrate  the  ability  to  generate  a  range  of  readily 
available  sets  of  constraint  orthogonal  polynomial  basis  functions. 


fix  =  .25)  =  o, 
fix  =  .25)  =  0 


0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1 


X - ► 


X - ► 


(a)  (b) 

Figure  3.  Constraint  Polynomial  Basis  Functions,  (a)  Basis  Function  are  Constraint  to  be  Zero  atx  = 
0.25.  (b)  Both,  Function  and  its  Derivative  are  Constrained  to  be  Zeros  atx  =  0.25 


B.  Ortho-Normal  Polynomials  Basis  Functions  in  Two-Dimensions 
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In  this  section  we  demonstrate  an  approach  to  produce  a  set  of  polynomial  orthonormal  basis  functions 
in  two  dimensions,  represented  by  X  and  y  with  maximum  degree  of  N,M  respectively.  First  we 
introduce  the  vector  space  of  two-dimensional  polynomials  and  then  the  algorithm  to  generate  orthonormal 
basis  vectors. 

1.  Matrix  Representation  of  2-D  Polynomials 


Let  PN>M  be  a  set  of  polynomials  where  P{x,y)^VN  M  is  a  polynomial  in  x  and  y  of  degree  N  and  M 


respectively.  Let  the  coefficients  be  a.. . 


P(x>y)-aNMxN yM  +aN_lMxN  lyM  +  ---  +  a0Mx°yM  +aN_lMxN  1  yM  +  a00x° y°  , 


or 


vA r-i  ^ M-j 

^7=1  W 

Let  us  represent  the  polynomial  coefficients  in  matrix  form  as, 


(13) 


aNM  aN-\M  aN-2M  ' ' '  a0 M 
aNM-\  aN-\M-\  aN-2M-\  a0M-\ 


(14) 


aN  0  aN- 10  aN- 2,0 


^00 


\m+i,n+i 


Then  the  polynomial  P(x,y)  can  be  represented  in  matrix  form  as, 

P{x,y)  =  XT  AY , 

where, 


~xN  1 

~yM~ 

x  = 

xN~l 

,  and 

Y  = 

yM- 1 

1 

N+\ 

1 

(15) 


(16) 


2.  Representation  of  a  2-D  Polynomial  Inner  Product  Vector  Space 

Let  PN  M  be  a  set  of  polynomials  in  x  and  y  of  degrees  up  to  N  and  M  respectively.  The  polynomial  set 
P N)m  defines  a  polynomial  vector  space  over  R  (set  of  real  numbers)  with 

®  =  (17) 

as  the  basis  vectors  that  satisfy  the  basic  properties  of  associativity,  commutativity,  existence  of  identity 
and  inverse,  and  properties  of  scalar  multiplication.  In  order  to  define  orthogonality,  we  must  define  an 
inner  product  space.  An  inner  product  space  is  a  vector  space  with  an  additional  structure  of  inner  product. 
We  define  the  inner  product  on  the  2-D  polynomial  vector  space  as, 

i  i 

(P(x,  y )}  Q(x, y ))  =  |  J P(x,  y )  0  Q(x, y )  dx dy ,  P(x,y),Q(x,y)  e  VN>M .  (18) 

oo 

The  polynomial  function  multiplication  (8) ,  inside  the  double  integrals  can  be  implemented  in  MATLAB 
using  2D  convolution  of  the  corresponding  coefficient  matrices  (see  function  conv2). 
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Our  goal  is  to  construct  an  orthonormal  basis  for  the  polynomial  vector  space  PN>M  .  As  a  first  step  we 
compute  an  orthonormal  basis  set  of  vectors  using  the  Gram-Schmidt  procedure.  Then  we  compute  the 
coefficients  of  a  polynomial  function  P{x,y)  in  the  new  transformed  space  with  orthonormal  basis. 

3.  The  Gram-Schmidt  Procedure 


If  V  is  a  set  of  L  linearly  independent  vectors,  then  the  Gram  Schmidt  procedure  can  be  used  to  find  set 
W  of  L  orthonormal  vectors  that  span  the  same  vector  space  as  V.  In  the  case  of  the  2D  polynomial  vector 
space,  we  have  a  set  of  (N+l)x(M+l)  linearly  independent  vectors  0  =  | 0i>J\  that  span  the  space  of  NxM 

degree  polynomials.  We  can  apply  Gram-Schmidt  procedure  to  this  vector  space  to  compute  the  set  of 
(N+l)x(M+l)  orthonormal  vectors  O  =  [^  ]. 


where, 


0ij 


^oo  =  h 

<Poi  =xly°-(xly°>& oo ) 

<P02  =X2y°  ~  (x2y°  ,000  )  ~  (X  V  >  0O\ )> 


(Pij=xJyl -(xJy\0j_u)-(xJylf0j_2^ - (x-y  Ao) 

and,  \Vij\  =  <\(<Pij><Pij)- 

Once  we  have  the  transformation  from  0  to  O  ,  we  can  transform  any  given  vector  (polynomial  in  this 
case)  to  the  new  orthonormal  basis  and  compute  the  new  set  of  coefficients,  by  taking  the  projection  of  the 
given  vector  (polynomial)  onto  each  of  the  new  basis  vectors. 

The  polynomial  vector  P(x,  j/)can  be  represented  in  a  new  basis  with  coefficients  • },  where 

aij  =  projection{P( x,y)^j)=<  P( x.y),^  >  . 

C.  Constrained  Orthonormal  Basis  Set  for  2-D  Functions 

Many  times  we  would  like  to  represent  a  class  of  functions  in  terms  of  basis  functions  that  must  satisfy 
certain  constraints.  For  an  efficient  representation,  we  would  like  to  construct  a  set  of  basis  functions 
satisfying  the  given  constraints  in  addition  to  the  orthonormality  condition.  In  this  report  we  restrict  our 
discussion  to  three  types  of  constraints;  constraints  on  the  function  value ,  its  first  derivative  or  both.  In  the 
following  we  describe  the  generation  of  2-D  constrained  basis  functions. 

(i)  Constraint  only  on  the  value  of  the  function  at  a  given  point.  I.e.,  C(x,y)  =  0  .  To  introduce 
the  constraints,  first  we  modify  the  polynomial  basis  set  0  by  multiplying  each  basis  function 
Oi  j{x,y) by  C{x,y )  to  form  a  new  constrained  basis  set  0C  given  as 

®c  =  \0ij(x>y)®c(x>y)}- 

Note  that  the  basis  set  0C  is  not  yet  orthonormal.  The  constrained  orthonormal  basis  set  Oc  is 
constructed  using  the  Grahm-  Schmidt  procedure  as  discussed  in  the  previous  section.  See  Figure 
4  for  an  example  of  the  orthonormal  basis  functions  generated. 

(ii)  Constraint  on  the  value  and  derivative  along  a  curve  C(x0,y0)  =  0  .  I.e.  C(x0,y0)  =  0  and 
C( x0,y0J  =  0  .  C  is  the  derivative  in  the  direction  orthogonal  to  the  curve  C(x0,y0)  =  0  in  the 
x-y  plane.  To  introduce  the  constraints,  we  first  define  C{x,y)2  =C(x,y)®C{x,y) .  Next  we 
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modify  the  polynomial  basis  set  0  by  multiplying  each  basis  function  Oij( x,y)  by  C( x,y )2  to 
form  a  new  constrained  basis  set  0C  given  as 

®cJpyx,y)®C(x,yf\. 

The  constrained  orthonormal  basis  set  Oc  is  constructed  using  the  Grahm-  Schmidt  procedure  as 
discussed  in  the  previous  section.  See  Figure  5  for  an  example  of  the  orthonormal  functions 
generated. 

(iii)  Constraint  only  on  the  derivative  of  the  function.  We  then  modify  the  given  polynomial  0  as, 
©c  =  Ko(x^)}^{©-Ko(x^)}}®C(x,j)2 . 

Figure  6  gives  an  example  of  orthonormal  functions  generated  with  constraints  only  on  the  derivatives. 


Figure  4.  Orthonormal  basis  function  set  of  polynomials  in  x  and  y  of  degrees  2  and  2,  with  the 
constraint  that  at  x+y  =  1,  the  value  of  the  function  is  zero. 
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Figure  5.  Orthonormal  basis  function  set  of  polynomials  in  x  and  y  of  degrees  2  and  2,  with  the 
constraint  that  at  x+y  =  1,  the  value  and  derivatives  of  the  function  are  zero. 


Figure  6.  Orthonormal  basis  function  set  of  polynomials  in  x  and  y  of  degrees  2  and  2,  with  the 
constraint  that  at  x+y  =  1,  the  derivative  of  the  function  is  zero. 
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VI.  Aerodynamic  Coefficient  Function  Estimation 

An  important  goal  of  the  ATRC  system  is  the  adaptive  reshaping  of  the  RLV  trajectory  in  the  presence 
of  altered  dynamic  characteristics  of  the  vehicle  when  unexpected  damage  occurs  in  the  various  operating 
scenarios.  Any  damage  to  a  vehicle  that  has  an  impact  on  the  external  shape  of  the  vehicle,  or  that  creates 
an  impediment  in  normal  functioning  of  the  control  surfaces,  results  in  alteration  of  the  vehicle’s 
aerodynamic  characteristics.  Figure  9  and  Figure  10  show  few  examples  of  the  Pitch  moment  coefficient 
variation  in  the  presence  of  various  damage  scenarios.  Since  the  aerodynamic  behavior  of  the  vehicle  is 
captured  in  aerodynamic  coefficients  that  are  used  for  the  design  of  vehicle  control  and  trajectory  planning, 
it  becomes  mission  critical  to  adapt  reference  trajectory  for  the  altered  vehicle  dynamics.  Hence  we  need  an 
approach  to  build  a  smooth  on-line  aerodynamic  model.  Physics  based  modeling  such  as  DATCOM  allows 
fast  computation  of  aerodynamic  coefficients  for  given  flight  points  and  the  results  can  be  stored  in  a 
tabular  form.  However,  for  efficient  real-time  trajectory  reshaping  application,  it  is  desired  to  represent 
aerodynamic  coefficients  in  smooth  functional  form  that  are  governed  by  few  parameters.  In  this  section, 
we  present  a  piecewise  continuous  and  smooth  function  model  in  two  dimensions. 

A.  Physic  Based  Modeling  for  Aerodynamic  Modeling 

Given  the  geometry  of  the  vehicle,  very  good  estimations  of  the  aerodynamic  coefficients  can  be 
generated  based  on  Physics  based  modeling  in  a  very  small  time.  For  example,  if  damage  on  the  vehicle  is 
characterized  adequately,  the  DATCOM  technology  allows  specification  of  the  altered  geometry  of  the 
vehicle  to  compute  the  corresponding  aerodynamic  coefficients  at  desired  points  of  the  flight  envelope.  In 
our  approach,  we  generate  a  large  number  of  data  points  for  the  aerodynamic  coefficients  and  then  fit  a 
piecewise  continuous  and  smooth  function  to  create  a  functional  model  for  vehicle  aerodynamics.  The 
functional  form  is  later  used  for  trajectory  reshaping. 

B.  Finite  Element  Function  Approximation 

We  formulate  the  aerodynamic  coefficient  function  with  a  set  of  parameters.  First  we  show  our 
approach  for  piecewise  continuous  and  smooth  one-dimensional  function,  which  is  later  extended  to  two- 
dimensions.  We  used  a  finite  element  modeling  approach  so  that  the  approximation  function  would  capture 
local  variations  in  an  efficient  manner.  Ref.  20,21  demonstrates  the  use  of  finite  element  piecewise 
approximation  for  mapping  geopotential.  Ref.  4,7  applied  the  technique  for  aerodynamic  coefficients 
representation.  First,  the  argument  space  of  the  function  is  divided  to  form  a  grid  with  one  control  point  at 
every  grid  point.  At  each  control  point  we  use  a  local  polynomial  function  that  is  determined  using  a 
weighted,  least  square  method  from  a  given  set  of  nominal  data  generated  from  DATCOM. 

4.  One-Dimensional  Finite  Element  Approximation 

The  scope  of  each  local  polynomial  centered  at  the  control  point  lies  in  between  the  adjacent  grid  points 
(see  Figure  7).  Notice  the  overlapping  of  local  polynomials,  which  helps  in  obtaining  a  smooth  function 
over  the  entire  range.  Once  local  approximations  are  determined,  a  smooth  global  approximation  is 
obtained  as  a  weighted  combination  of  the  local  approximations.  The  smoothness  of  the  approximation 
function  implies  that  the  function,  as  well  as  its  first  derivative,  is  continuous.  Notice  that  grids  are  not 
restricted  to  be  equidistant.  To  capture  nonlinearity  effectively,  more  control  points  should  be  placed  near 
high  non-linearity.  Figure  7  uses  a  second  order  polynomial  function  for  local  approximations. 

Once  local  approximations  are  attained,  a  smooth  global  approximation  is  obtained  as  a  weighted 
combination  of  the  local  approximations.  A  highlight  of  this  Finite  Element  approximation  is  that  it 
preserves  the  local  function  value  and  its  first  derivative  at  its  control  point.  This  is  achieved  by  a  smooth 
weighting  function  that  gracefully  goes  from  unity  to  zero,  from  one  control  point  to  another,  without 
contributing  to  a  first  derivative  at  both  control  points.  The  weighting  function  for  the  two  overlapping 
local  approximation  curves  is  given  as 

W,(q)  =  {\-qf{\  +  2q),  W2(q)  =  \-Wl(q)  , 

where  q  is  the  normalized  coordinate  with  the  origin  at  the  first  control  point.  Note  that  the  weighting 
function  ensures  that  Wl(0)  =  l9  W2(0)  =  0  at  the  first  control  point  and  similarly,  W1(1)  =  0,  W2(l)  =  l 
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at  the  second  control  point.  If  f\(q)  and  /2  (q )  are  two  local  approximations,  the  weighted  global 
approximation  F(q)  between  two  control  points  is  given  as 

F(q)  =  W,  (q )fi (q)  +  W2 (q)f2 (q)  . 

Figure  8  shows  the  final  approximated  function  that  is  made  of  the  weighted  combinations  of  the  local 
function  approximations. 


Figure  7:  Local  Approximations  Centered  at  Figure  8:  Smooth  Functional  Approximation 

Control  Points  using  Weighted  Local  Approximations 
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Figure  9.  Moment  Coefficient  with  AOA  for 
Nominal  and  Various  Failure  Cases 


Figure  10.  Moment  Coefficient  with  Mach  for 
Nominal  and  Various  Failure  Cases 


5.  Multi-Dimensional  Finite  Element  Approximation 
In  this  section  we  extend  the  formulation  for  FEM  modeling  for  multi  dimensional  functions.  Let  F  be 
given  as 

F  =  F(ql,q2,  —  ,qH)  *  (19) 

and  we  have  to  determine  an  estimate  F(ql,q2p  •  -,qn)  from  a  finite  element  model.  As  in  one¬ 
dimensional  case,  we  assume  that  the  domain  of  F  is  covered  with  finite  number  of  equidistant  nodes. 
Each  preliminary  local  approximation  Ft  ^  .  is  valid  in  the  2  x  2  x  •  •  •  x  2  hypercube  centered  at  a  node 

represented  by  n  indices  given  as  ( z t ,  z2 ,  ■  •  * ,  )  .  Each  index  is  numbered  in  increasing  order  for  the 

nodes  lying  in  its  dimension.  Consider  a  hyper  cube  formed  by  2n  nodes,  each  node  representing  a  comer 
of  the  hyper  cube. The  final  approximation  Ft  ^  .  is  valid  in  the  unit  hypercube  whose  “lower  left 
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comer”  is  (/, ,  i2,  •  •  ■ , in  )  .  We  shift  the  origin  of  the  coordinates  to  the  node  (7, ,  /2 ,  •  •  • , in )  that  has  the 
lowest  indices  in  each  dimension.  The  local  coordinates  are  normalized  as 


cli  =  - 


Vj-VU: 


y'd'y+1 


Let  us  define  a  function  W  as  a  function  of  the  normalized  coordinate  qj  as 

W  (q, )  =  (1  -  iL  )2  (1  +  2q. ) 


(20) 


(21) 


Now  the  weight  function  for  an  arbitrary  comer  node  (7, ,  /2 ,  •  •  • ,  ik  + 1,  ik+l  + 1,  •  •  • ,  in  )  of  the  hyper  cube 


can  be  written  as 

"U  ,,4^,+v.v,  =  )W(q2  )-W(l-qt  )W(  1  -  ft.,)  •  •  •  ? 

for  0  <  q.  <  1,  j  =  l-n  (22) 

The  weight  function  defined  by  Eq.  (  22)has  the  following  properties: 

•  weight  is  unity  at  its  own  node  and,  it  is  zero  at  any  other  node  or  edge  formed  by  other 
nodes. 

•  summation  of  all  the  weights  corresponding  to  the  comer  nodes  of  the  hyper  cube  at  any 
point  inside  the  hyper  cube  is  equal  to  one. 

The  estimated  function  F  can  be  defined  in  terms  of  normalized  coordinates  centered  at  node 


(  h  ?  4  ?  ‘  ?  4  ) 


as 


F(ql,q2,---,q„)  =  'Yj 


w. 


■>W  *1»*2  >-»W 


j=1  .  (23) 

6.  An  Example  of  Finite  Element  Function  Approximation  for  Functions  of  Two  Variables. 

Let  Z  be  the  matrix  of  values  of  a  function  sampled  at  the  grid  points  of  a  2-d  space  defined  by  x  andy.  We 
are  interested  in  representing  the  function  Z  in  closed  form,  using  2-d  polynomial  approximation.  The 
closed  form  representation  can  be  used  to  evaluate  the  function  Z  at  arbitrary  x  and  y  points.  We 
approximate  the  function  Z  using  a  number  of  2-d  polynomial  function  elements  of  a  smaller  area  of 
support.  The  number  of  function  elements  to  be  used  and  the  degrees  of  the  polynomials  in  x  and  y  are 
inputs  to  the  algorithm.  We  specify  this  by  the  x  andy  coordinates  of  the  control  points  as  separate  vectors. 
Note  that  the  control  points  lie  exactly  on  the  grid  points.  At  these  data  points,  the  approximated  value  of 
the  function  Z  equals  the  value  of  the  local  function  centered  at  that  control  point.  The  chosen 
approximation  is  locally  optimal  in  the  least  squared  error  criteria  (i.e.,  within  the  area  of  support  of  each 
function  element). 

To  formulate  finite  element  function  approximation  satisfying  certain  constraints,  we  proceed  as  follows. 
The  orthonormal  basis  set  of  functions  are  computed  that  satisfy  the  constraints  and  then  searched  for  local 
functions  that  are  optimal  in  the  vector  space  spanned  by  the  constrained  orthonormal  basis  functions, 
based  on  least  square  optimality  criteria.  More  details  on  the  finite  element  function  approximation  can  be 
found  in  [7]. 


Concentrating  on  the  finite  element  function  approximation  method  for  2-d  functions,  we  initially 
proceeded  by  calculating  the  weights  in  the  x  and  y  directions  separately.  The  actual  weight  function  at 
each  point  in  the  grid  is  the  product  of  the  x  and  y  weights  at  that  location.  The  x  and  y  weight  functions 
satisfy  the  criteria  that  at  the  current  control  point  location,  it  has  a  value  of  one  and  gradually  decreases  to 
zero  as  it  moves  away.  The  products  of  the  1-d  weight  functions  also  satisfy  the  same  criteria  in  the  2-d  x-y 
plane.  The  intended  consequence  of  this  property  is  that  the  approximated  function  value  will  be  equal  to 
the  value  of  the  local  function  at  the  control  point.  Figure  1 1  shows  an  example  of  a  weight  function  in  the 
2-d  plane  that  is  the  product  of  the  two  1-d  weight  functions. 
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Figure  11.  Weight  function  centered  at  a  control  point.  The  value  is  one  at  the  location  of  the 
control  point  and  gradually  decreases  and  reaches  zero  at  the  neighboring  control  points. 


Figure  12  shows  the  results  of  finite  element  approximation  for  a  surface  of  31x31  points  with  5x5  =  25 
control  points  and  local  polynomials  of  degree  3  in  x  and  y. 


Z-  Given  function  FEM  approximation  ofZ  Error  in  function  approximation 


Figure  12.  Finite  element  function  approximation,  (a)  Given  function  (b)  Approximated  Function  (c) 
Error  between  the  given  and  approximated  functions.  Twenty-five  control  points  were  used  and 
each  local  function  is  a  polynomial  of  degree  3  in  x  and  3  in  y. 
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VII.  Conclusion 

In  this  paper  we  address  the  function  modeling  techniques  to  facilitate  an  efficient  on-line  inverse 
dynamics  methodology  for  trajectory  reshaping  of  the  RLVs  .  On-line  trajectory  reshaping  to  determine  a 
feasible  reference  trajectory  is  computationally  a  difficult  problem  for  real  time  applications. 
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